package com.radar.task;

import cn.hutool.core.date.DatePattern;
import cn.hutool.core.date.DateTime;
import cn.hutool.core.date.DateUtil;
import cn.hutool.core.io.FileUtil;
import cn.hutool.core.util.ZipUtil;
import com.google.gson.Gson;
import com.radar.algorithm.*;
import com.radar.common.*;
import com.radar.controller.RadarController;
import com.radar.exception.RadarCaclError;
import com.radar.exception.RadarDecodeError;
import com.radar.io.CinradRadar;
import com.radar.io.RadarGridData;
import com.radar.io.cinrad.RadarCC;
import com.radar.utils.BitConverter;
import com.radar.utils.RefObj;
import com.radar.visualize.PolarMap;
import com.radar.visualize.RasterMap;
import org.apache.commons.io.filefilter.RegexFileFilter;
import org.bytedeco.javacpp.indexer.Indexer;
import org.bytedeco.opencv.global.opencv_imgcodecs;
import org.bytedeco.opencv.global.opencv_imgproc;
import org.bytedeco.opencv.opencv_core.Mat;
import org.locationtech.jts.geom.Coordinate;
import org.meteoinfo.global.Extent;
import org.meteoinfo.global.PointD;
import org.springframework.beans.factory.annotation.Value;
import org.springframework.scheduling.annotation.Scheduled;
import org.springframework.stereotype.Component;
import ucar.ma2.InvalidRangeException;

import javax.imageio.ImageIO;
import java.awt.*;
import java.awt.geom.Line2D;
import java.awt.image.BufferedImage;
import java.io.File;
import java.io.IOException;
import java.nio.charset.Charset;
import java.util.*;
import java.util.List;

@Component
public class RocketForceTask extends BaskTask {

    @Value("${radar.root-path:''}")
    public void setRootPath(String rootPath) {
        super.rootPath = rootPath;
    }

    @Value("${radar.process-path:''}")
    public void setProcessPath(String processPath) {
        super.processPath = processPath;
    }

    @Value("${radar.handle-hour:-1}")
    public void setHandleHour(int handleHour) {
        super.handleHour = handleHour;
    }

    @Value("${radar.radar-station:''}")
    private String radarcode;

    @Override
//    @Scheduled(cron = "0 0/1 * * * ? ")
    @Scheduled(fixedDelay = 1000 * 60)
    public void jobRun() {
        if (!super.filterDir()) {
            return;
        }
        log("start radar task poll date : " + DateUtil.formatDateTime(DateUtil.date()), "RADAR-SINGLE");
        long time = System.currentTimeMillis();
        FileUtil.mkdir(processPath);
        FileUtil.loopFiles(this.rootPath, new RegexFileFilter("^Z_RADR_I_Z(\\d{4})_(\\d{14})_O_DOR_(.*?)_CAP_FMT.(BIN.BZ2|bin.bz2|bin|BIN)$"))
                .stream()
                .sorted(Comparator.comparingLong(i -> i.lastModified()))
                .filter(i -> {
                    if (this.handleHour == -1) {
                        return true;
                    } else {
                        // handle hour
                        Date date = DateUtil.parse(i.getName().split("_")[4].substring(0, 12), DateConstants.DATE_FORMAT_YYYYMMDDHHMM);
                        return date.getTime() > DateUtil.offsetHour(new Date(), -handleHour).getTime();
                    }
                }).forEach(i -> {
                    try {
                        this.process(i);
                    }catch (Exception e){
                        e.printStackTrace();
                    }
                });
        log("end radar task poll times: " + (System.currentTimeMillis() - time), "RADAR-SINGLE");
    }

    @Override
    public void process(File file) {
        File OKFile = new File(super.okPath + File.separator + "RADAR-SINGLE" + File.separator + file.getName() + ".OK");
        if (OKFile.exists()) {
            return;
        } else {
            FileUtil.touch(OKFile);
        }
        log("start parse radar file : " + file.getAbsolutePath(), "RADAR-SINGLE");
        DateTime parse = DateUtil.parse(file.getName().split("_")[4], DatePattern.PURE_DATETIME_PATTERN);
        String relativePath = this.processPath + File.separator + DateUtil.format(parse, DateConstants.DATE_FORMAT_YYYYMMDD) + File.separator + DateUtil.format(parse, DateConstants.DATE_FORMAT_YYYYMMDDHHMM) + File.separator;
        FileUtil.mkdir(relativePath);
        // cp file
        File target = new File(relativePath + file.getName());
        FileUtil.copy(file, target, true);
        try {
            long time = System.currentTimeMillis();

            CinradRadar cinradRadar = new CinradRadar();
            // TODO 修改自定义站点
            if (radarcode.equals("")) {
                cinradRadar.parse(target.getAbsolutePath());
            } else {
                cinradRadar.parse(target.getAbsolutePath(), null, radarcode);
            }
            log("start parse radar file " + target.getName(), "RADAR-SINGLE");
            cinradRadar.radar.setProjectStep(0.02);

            Date date = cinradRadar.radar.dateTime;
            RadarSite.RadarSiteItem radarSite = cinradRadar.radar.getRadarSite();
            Map<String, Object> result = new HashMap<>();
            result.put("code", cinradRadar.code);
            result.put("alt", cinradRadar.alt);
            result.put("radarSite", radarSite);
            result.put("date", DateUtil.format(cinradRadar.radar.dateTime, DateConstants.DATE_TIME_FORMAT_YYYY_MM_DD_HH_MI_SS));
            result.put("type", cinradRadar.type.name());
            result.put("dataType", cinradRadar.datType.name());
            result.put("vcp", cinradRadar.radar.vcp);
            result.put("elevation", cinradRadar.radar.getElevation());
            result.put("rresolution", cinradRadar.radar.getResolution(ProductType.REF));
            result.put("vresolution", cinradRadar.radar.getResolution(ProductType.VEL));
            // 写入基本信息
            FileUtil.writeUtf8String(new Gson().toJson(result), relativePath + "setting.json");
            // 画图
            int[] refIndex = cinradRadar.radar.getElevationIndex(ProductType.REF);
            int[] velIndex = cinradRadar.radar.getElevationIndex(ProductType.VEL);
            BufferedImage bufferedImage;

            log("开始处理 单层 ppi station: " + radarSite.getStation() + ",date : " + DateUtil.format(date, DateConstants.DATE_FORMAT_YYYYMMDDHHMM) + ",path : " + target.getAbsolutePath(), "RADAR-SINGLE");
            ImageIO.setUseCache(false);
            PolarMap polarMap = new PolarMap(cinradRadar.radar);
            RasterMap rasterMap = new RasterMap(cinradRadar.radar);
            for (int i = 0; i < refIndex.length; i++) {
//                bufferedImage = rasterMap.simplePPI(RasterMap.buildRasterInfo(ProductType.REF)
//                        .simpleDetail(true)
//                        .drawCity(true), refIndex[i]);
                bufferedImage = rasterMap.simplePPI(RasterMap.buildRasterInfo(ProductType.REF).backColor(new Color(0, 0, 0, 0)), refIndex[i]);
                ImageIO.write(bufferedImage, "png", new File(relativePath + "REF_" + i + ".png"));

                bufferedImage = rasterMap.simplePPI(RasterMap.buildRasterInfo(ProductType.VEL).backColor(new Color(0, 0, 0, 0)), velIndex[i]);
                ImageIO.write(bufferedImage, "png", new File(relativePath + "VEL_" + i + ".png"));
            }
            // draw cr
            log("开始处理 et ppi station: " + radarSite.getStation() + ",date : " + DateUtil.format(date, DateConstants.DATE_FORMAT_YYYYMMDDHHMM) + ",path : " + target.getAbsolutePath(), "RADAR-SINGLE");

            RadarGridData et = cinradRadar.radar.getET(0);
            bufferedImage = rasterMap.simpleRaster(RasterMap.buildRasterInfo(ProductType.ET).backColor(new Color(0, 0, 0, 0)), et);
            ImageIO.write(bufferedImage, "png", new File(relativePath + "ET.png"));

            log("开始处理 vil ppi station: " + radarSite.getStation() + ",date : " + DateUtil.format(date, DateConstants.DATE_FORMAT_YYYYMMDDHHMM) + ",path : " + target.getAbsolutePath(), "RADAR-SINGLE");
            RadarGridData vil = cinradRadar.radar.getVIL(0);
            bufferedImage = rasterMap.simpleRaster(RasterMap.buildRasterInfo(ProductType.VIL).backColor(new Color(0, 0, 0, 0)), vil);
            ImageIO.write(bufferedImage, "png", new File(relativePath + "VIL.png"));

            log("开始处理 ppi_ref ppi station: " + radarSite.getStation() + ",date : " + DateUtil.format(date, DateConstants.DATE_FORMAT_YYYYMMDDHHMM) + ",path : " + target.getAbsolutePath(), "RADAR-SINGLE");
            bufferedImage = rasterMap.simplePPI(RasterMap.buildRasterInfo(ProductType.REF).backColor(new Color(0, 0, 0, 0)), cinradRadar.radar.getElevationIndex(ProductType.REF)[0]);
            ImageIO.write(bufferedImage, "png", new File(relativePath + "PPI_REF.png"));

            log("开始处理 ppi_vel ppi station: " + radarSite.getStation() + ",date : " + DateUtil.format(date, DateConstants.DATE_FORMAT_YYYYMMDDHHMM) + ",path : " + target.getAbsolutePath(), "RADAR-SINGLE");
            bufferedImage = rasterMap.simplePPI(RasterMap.buildRasterInfo(ProductType.VEL).backColor(new Color(0, 0, 0, 0)), cinradRadar.radar.getElevationIndex(ProductType.VEL)[0]);
            ImageIO.write(bufferedImage, "png", new File(relativePath + "PPI_VEL.png"));


            log("开始处理 cr station: " + radarSite.getStation() + ",date : " + DateUtil.format(date, DateConstants.DATE_FORMAT_YYYYMMDDHHMM) + ",path : " + target.getAbsolutePath(), "RADAR-SINGLE");
            RadarGridData radardata = cinradRadar.radar.getCR();
            RasterMap.RasterInfo rasterInfo = RasterMap.buildRasterInfo(ProductType.CR).backColor(new Color(0, 0, 0, 0));
            bufferedImage = rasterMap.simpleRaster(rasterInfo, radardata);
            ImageIO.write(bufferedImage, "png", new File(relativePath + "CR.png"));

            // 计算冰雹
//            log("开始处理 hcl station: " + radarSite.getStation() + ",date : " + DateUtil.format(date, DateConstants.DATE_FORMAT_YYYYMMDDHHMM) + ",path : " + target.getAbsolutePath(), "RADAR-SINGLE");
//            cinradRadar.radar.getHCL();

            log("开始处理 gray station: " + radarSite.getStation() + ",date : " + DateUtil.format(date, DateConstants.DATE_FORMAT_YYYYMMDDHHMM) + ",path : " + target.getAbsolutePath(), "RADAR-SINGLE");
            double[][] gray = MathCommon.decode(radardata.radarData, ProductType.REF);
            gray = Cacl.nineMedianFilter(gray);
            bufferedImage = rasterMap.drawSingleRaster(gray, ProductType.REFGRAY);
            ImageIO.write(bufferedImage, "png", new File(relativePath + "GRAY.png"));

            // find forecast
            double[][][] dir = (double[][][]) null;
            String prePath = null;
            int offset = -1;
            for (int i = 1; i < 20; i++) {
                Date PreDate = DateUtil.offsetMinute(parse, -i);
                prePath = this.processPath + File.separator + DateUtil.format(PreDate, DateConstants.DATE_FORMAT_YYYYMMDD) + File.separator + DateUtil.format(PreDate, DateConstants.DATE_FORMAT_YYYYMMDDHHMM) + File.separator + "GRAY.png";
                if (FileUtil.exist(prePath)) {
                    log("start forecast prePath step : " + i, "RADAR-SINGLE");
                    dir = forecast(prePath, relativePath + "GRAY.png", i, radardata, relativePath);
                    offset = i;
                    break;
                }
            }
            // output storms 计算
            RefObj<Integer> count1 = new RefObj<>(0);
            int step = 8;
            double[][] decode = MathCommon.decode(radardata.radarData, ProductType.REF);
            double[][] data = new double[decode.length / step][decode[0].length / step];
            // 降低分辨率 不然算的太慢
            for (int i = 0; i < data.length; i++) {
                for (int j = 0; j < data[0].length; j++) {
                    data[i][j] = decode[i * step][j * step];
                }
            }

            step = 4;
            data = new double[decode.length / step][decode[0].length / step];
            for (int i = 0; i < data.length; i++) {
                for (int j = 0; j < data[0].length; j++) {
                    data[i][j] = decode[i * step][j * step];
                }
            }
            List<Storms.StormDt> list = Storms.getStorms(data, step, count1, 32, 1000, cinradRadar);

            Gson gson = new Gson();
            double[][][] finalDir = dir;
            // 用光流替代titan
            if (finalDir != null) {
                // 写入 uv 图片
                int h = Math.min(finalDir[0].length, decode.length);
                int w = Math.min(finalDir[0][0].length, decode[0].length);
                // 生成风场
                bufferedImage = rasterMap.drawSingleRaster(MathCommon.decode(radardata.radarData, ProductType.REF), ProductType.REF, new Color(0, 0, 0, 0));
                BufferedImage br = new BufferedImage(w, h, 2);
                Graphics2D graphics = (Graphics2D) bufferedImage.getGraphics();
                Graphics2D graphics1 = (Graphics2D) br.getGraphics();
                graphics.setColor(Color.black);
                graphics1.setColor(Color.black);
                for (int ii = 0; ii < h; ii += 25) {
                    for (int jj = 0; jj < w; jj += 25) {
                        double x = finalDir[0][ii][jj];
                        double y = finalDir[1][ii][jj];
//                    if (decode[ii][jj] > 5) {
                        double ws = Math.sqrt(x * x + y * y);
                        WindFlagSymbol wfs = new WindFlagSymbol(ws);
                        wfs.rotateTo(x, y);
                        List<Line2D> awtLines = wfs.getAwtLines(new Coordinate(jj, ii), 15);
                        awtLines.forEach(graphics::draw);
                        awtLines.forEach(graphics1::draw);

                        graphics.setFont(new Font("Arial", 0, 10));
                        graphics1.setFont(new Font("Arial", 0, 10));
                        graphics.drawString((int) WindFlagSymbol.getDirAndSpeed(x, y) + "", jj, ii);
                        graphics1.drawString((int) WindFlagSymbol.getDirAndSpeed(x, y) + "", jj, ii);
                    }
                }
                ImageIO.write(bufferedImage, "png", new File(relativePath + File.separator + "CR_UV.png"));
                ImageIO.write(br, "png", new File(relativePath + File.separator + "CR_UV_BLANK.png"));

                int[][][] d = new int[finalDir.length][finalDir[0].length][finalDir[0][0].length];
                for (int ii = 0; ii < h; ii ++) {
                    for (int jj = 0; jj < w; jj ++) {
                        if (decode[ii][jj] > 5) {
                            double x = finalDir[0][ii][jj];
                            double y = finalDir[1][ii][jj];
                            double ws = Math.sqrt(x * x + y * y);
                            d[0][ii][jj] = (int) (ws * 10);
                            d[1][ii][jj] = (int) WindFlagSymbol.getDirAndSpeed(x, y);
                        }
                    }
                }
                String uvpath = relativePath + File.separator + "uv.dat";
                File f = new File(uvpath);
                FileUtil.writeBytes(BitConverter.getBytes(d), f);
                ZipUtil.zip(f.getAbsolutePath(), f.getAbsolutePath() + ".radar");
                f.delete();

                Map<String, Object> uvresult = new HashMap<>();
                uvresult.put("h", d[0].length);
                uvresult.put("w", d[0][0].length);
                // 写入基本信息
                FileUtil.writeUtf8String(new Gson().toJson(uvresult), relativePath + "uv_setting.json");

            }
            for (Storms.StormDt z : list) {
                if (finalDir != null) {
                    double u = 0.0D;
                    double v = 0.0D;
                    int xindex = (int) (z.getCenter()).Y;
                    int yindex = (int) (z.getCenter()).X;
                    int index = 0;
                    for (int i = yindex - 4; i <= yindex + 4; i++) {
                        for (int j = xindex - 4; j <= xindex + 4; j++) {
                            if (decode[i][j] > 5 && j > 0 && i > 0 && i < finalDir[0].length - 1 && j < finalDir[0][0].length - 1) {
                                u += finalDir[0][i][j];
                                v += finalDir[1][i][j];
                                index++;
                            }
                        }
                    }
                    if (index > 0) {
                        u /= index;
                        v /= index;
                        z.setWd(WindFlagSymbol.getDirAndSpeed(u, v));
                        z.setWs(cinradRadar.radar.getResolution(ProductType.CR) / offset / 60 * uv2wind(u, v)[0]);
                    }
                }
                if (z.getWs() != -999 && z.getWd() != -999) {
                    z.setForecastCenter(computerThatLonLat(z.getCenter(), z.getWd(), z.getWs() * 120 * 60));
                }
                List<double[]> list1 = new ArrayList<>();
                z.getBoundary().forEach((point) -> {
                    double x = radardata.xArray[(int) point.Y];
                    double y = radardata.yArray[radardata.yArray.length - 1 - (int) point.X];
                    point.X = x;
                    point.Y = y;
                    list1.add(new double[]{y, x});
                });
                z.setBoundary707(list1);
                z.getContain().forEach((point) -> {
                    double x = radardata.xArray[(int) point.Y];
                    double y = radardata.yArray[radardata.yArray.length - 1 - (int) point.X];
                    point.X = x;
                    point.Y = y;
                });
            }
            FileUtil.writeString(gson.toJson(list), new File(relativePath + "storms.json"), Charset.defaultCharset());
            FileUtil.writeString(radardata.getXMin() + ";" + radardata.getXMax() + ";" + radardata.getYMin() + ";" + radardata.getYMax(),
                    new File(relativePath + File.separator + "extent.txt"), Charset.defaultCharset());
            FileUtil.touch(OKFile);

            // write 3d
            String rPath = relativePath + "surface.dat";
            double[] alt = new double[]{1000, 1500, 2000, 2500, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000, 11000, 12000};
            byte[][][] cappi = cinradRadar.radar.cappi(ProductType.REF, alt);
            double jself = 500.0 / cappi[0].length;
            double kself = 500.0 / cappi[0][0].length;
            byte[] bytes = new byte[alt.length * 500 * 500];
            int index = 0;
            for (int i = 0; i < alt.length; i++) {
                bufferedImage = rasterMap.drawSingleRaster(MathCommon.decode(cappi[i], ProductType.REF), new Color(0, 0, 0, 0));
                ImageIO.write(bufferedImage, "png", new File(relativePath + "CAPPI_" + (int) alt[i] + ".png"));

                for (int j = 0; j < 500; j++) {
                    for (int k = 0; k < 500; k++) {
                        int jj = (int) Math.floor((500 - 1 - j) / jself);
                        int kk = (int) Math.floor(k / kself);
                        if (cappi[i][jj][kk] == -128 || MathCommon.DecodeRef(cappi[i][jj][kk]) <= 5) {
                            bytes[index] = (byte) 255;
                        } else {
                            bytes[index] = (byte) MathCommon.DecodeRef(cappi[i][jj][kk]);
                        }
                        index++;
                    }
                }
            }
            File f = new File(rPath);
            FileUtil.writeBytes(bytes, f);
            ZipUtil.zip(f.getAbsolutePath(), f.getAbsolutePath() + ".zip");
            f.delete();
            cinradRadar.radar.flush();
            log("end parse radar file " + target.getName() + " time : " + (System.currentTimeMillis() - time), "RADAR-SINGLE");
        } catch (Exception e) {
            log("error write radar file : " + file.getName() + " e : " + e.getMessage(), "RADAR-SINGLE");
        }
    }

    public Extent getExtent(ProductType type, CinradRadar cinradRadar, int index) throws RadarCaclError, RadarDecodeError, IOException, InvalidRangeException {
        Extent extent = null;
        double lon = cinradRadar.radar.getRadarSite().getLon();
        double lat = cinradRadar.radar.getRadarSite().getLat();
        if (type == ProductType.REF || type == ProductType.VEL) {
            double res = cinradRadar.radar.getResolution(type);
            byte[][] raw = cinradRadar.radar.getRaw(index, type);
            double ele = cinradRadar.radar.getElevation()[index];
            double distance = raw[0].length * res;
            distance = Math.cos(ele * Const.DEG2RAD) * distance;
            extent = Common.getExtentRadius(lon, lat, distance);
        } else if (type == ProductType.CR) {
            extent = cinradRadar.radar.getCR().getExtent();
        } else if (type == ProductType.ET) {
            extent = cinradRadar.radar.getET().getExtent();
        } else if (type == ProductType.EB) {
            extent = cinradRadar.radar.getEB().getExtent();
        } else if (type == ProductType.VIL) {
            extent = cinradRadar.radar.getVIL().getExtent();
        }
        return extent;
    }

    public double[][][] forecast(String preImg, String nextImg, int offset, RadarGridData cr, String output) throws IOException {
        if (!FileUtil.exist(preImg)) {
            return null;
        }
        Mat ff1 = opencv_imgcodecs.imread(preImg);
        Mat ff2 = opencv_imgcodecs.imread(nextImg);
        Mat prvs = ff1.clone();
        Mat next = ff1.clone();
        opencv_imgproc.cvtColor(ff1, prvs, opencv_imgproc.COLOR_BGR2GRAY);
        opencv_imgproc.cvtColor(ff2, next, opencv_imgproc.COLOR_BGR2GRAY);

        int step = 4;

        int w = Math.min(prvs.size().get(0), next.size().get(0)) / step;
        int h = Math.min(prvs.size().get(1), next.size().get(1)) / step;

        double[][] pre = new double[h][w];
        double[][] nex = new double[h][w];
        Indexer indexer = prvs.createIndexer();
        Indexer indexer1 = next.createIndexer();

        for (int i = 0; i < h; i++) {
            for (int j = 0; j < w; j++) {
                pre[i][j] = indexer.getDouble(i * step, j * step);
                if (pre[i][j] == 255) {
                    pre[i][j] = 0;
                }
                nex[i][j] = indexer1.getDouble(i * step, j * step);
                if (nex[i][j] == 255) {
                    nex[i][j] = 0;
                }
            }
        }


        double[][][] cacl = DenseLucaskanade.cacl(new double[][][]{pre, nex});
        double[][] uu = cacl[0];
        double[][] vv = cacl[1];
        for (int i = 0; i < h; i++) {
            for (int j = 0; j < w; j++) {
                uu[i][j] = uu[i][j] * (false ? 1.5 : 6.0 / offset);
                vv[i][j] = vv[i][j] * (false ? 1.5 : 6.0 / offset);
            }
        }

        uu = Cacl.nineMedianFilter(uu);
        uu = Cacl.nineMedianFilter(uu);
        vv = Cacl.nineMedianFilter(vv);
        vv = Cacl.nineMedianFilter(vv);

        double[][] u = new double[uu.length * step][uu[0].length * step];
        double[][] v = new double[uu.length * step][uu[0].length * step];
        for (int i = 0; i < u.length; i++) {
            for (int j = 0; j < u[0].length; j++) {
                u[i][j] = uu[i / step][j / step];
                v[i][j] = vv[i / step][j / step];
            }
        }

        double[][] decode = MathCommon.decode(cr.radarData, ProductType.REF);
        if (decode.length == u.length && decode[0].length == u[0].length) {

        } else {
            double[][] z = new double[u.length][u[0].length];
            for (int i = 0; i < u.length; i++) {
                for (int j = 0; j < u[0].length; j++) {
                    if (i >= decode.length || j >= decode[0].length) {
                        z[i][j] = GisConstants.UNDEF;
                    } else {
                        z[i][j] = decode[i][j];
                    }
                }
            }
            decode = z;
        }
        for (int i = 0; i < u.length; i++) {
            for (int j = 0; j < u[0].length; j++) {
                u[i][j] = uu[i / step][j / step];
                v[i][j] = vv[i / step][j / step];
            }
        }
        double[][][] dir = new double[2][u.length][(u[0]).length];
        for (int k = 0; k < u.length; k++) {
            for (int n = 0; n < (u[0]).length; n++) {
                dir[0][k][n] = u[k][n];
                dir[1][k][n] = v[k][n];
            }
        }
        SemiLagrangian semiLagrangian = new SemiLagrangian();
        int t = 20;
        float[][][] extrapolate = semiLagrangian.extrapolate(MathCommon.doubleToFloat(decode), MathCommon.doubleToFloat(u), MathCommon.doubleToFloat(v), t);

        for (int i = 0; i < t; i++) {
            BufferedImage bufferedImage;
            for (int j = 0; j < u.length; j++) {
                for (int k = 0; k < u[0].length; k++) {
                    extrapolate[i][j][k] = extrapolate[i][j][k] == 0 || extrapolate[i][j][k] == Float.NaN ? GisConstants.UNDEF : extrapolate[i][j][k];
                }
            }
            RasterMap rasterMap = new RasterMap(null);
            bufferedImage = rasterMap.drawSingleRaster(MathCommon.floatToDouble(extrapolate[i]), new Color(0, 0, 0, 0));
            ImageIO.write(bufferedImage, "png", new File(output + File.separator + "CR_" + (i + 1) + ".png"));
        }
        return dir;
    }

    public static double[] uv2wind(double u, double v) {
        double[] ws = new double[2];
        ws[0] = Math.sqrt(u * u + v * v);
        double dir = 270.0D - Math.atan2(v, u) * 180.0D / Math.PI;
        if (dir >= 360.0D)
            dir -= 360.0D;
        if (dir < 0.0D)
            dir += 360.0D;
        ws[1] = dir;
        return ws;
    }


    private static double a = 6378137;
    private static double b = 6356752.3142;
    private static double f = 1 / 298.2572236;

    public static PointD computerThatLonLat(PointD pointD, double brng, double dist) {
        double lon = pointD.X;
        double lat = pointD.Y;
        double alpha1 = rad(brng);
        double sinAlpha1 = Math.sin(alpha1);
        double cosAlpha1 = Math.cos(alpha1);
        double tanU1 = (1 - f) * Math.tan(rad(lat));
        double cosU1 = 1 / Math.sqrt((1 + tanU1 * tanU1));
        double sinU1 = tanU1 * cosU1;
        double sigma1 = Math.atan2(tanU1, cosAlpha1);
        double sinAlpha = cosU1 * sinAlpha1;
        double cosSqAlpha = 1 - sinAlpha * sinAlpha;
        double uSq = cosSqAlpha * (a * a - b * b) / (b * b);
        double A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
        double B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));

        double cos2SigmaM = 0;
        double sinSigma = 0;
        double cosSigma = 0;
        double sigma = dist / (b * A), sigmaP = 2 * Math.PI;
        while (Math.abs(sigma - sigmaP) > 1e-12) {
            cos2SigmaM = Math.cos(2 * sigma1 + sigma);
            sinSigma = Math.sin(sigma);
            cosSigma = Math.cos(sigma);
            double deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)
                    - B / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM)));
            sigmaP = sigma;
            sigma = dist / (b * A) + deltaSigma;
        }

        double tmp = sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1;
        double lat2 = Math.atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1,
                (1 - f) * Math.sqrt(sinAlpha * sinAlpha + tmp * tmp));
        double lambda = Math.atan2(sinSigma * sinAlpha1, cosU1 * cosSigma - sinU1 * sinSigma * cosAlpha1);
        double C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
        double L = lambda - (1 - C) * f * sinAlpha
                * (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));

        double revAz = Math.atan2(sinAlpha, -tmp); // final bearing
        return new PointD(lon + deg(L), deg(lat2));
    }

    private static double rad(double d) {
        return d * Math.PI / 180.0;
    }

    private static double deg(double x) {
        return x * 180 / Math.PI;
    }

}
